Shocks and excitation dynamics in a driven diffusive two-channel system 
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We consider classical hard-core particles hopping stochastically on two parallel chains in the same 
or opposite directions with an inter- and intra-chain interaction. We discuss general questions con- 
cerning elementary excitations in these systems, shocks and rarefaction waves. From microscopical 
1 considerations we derive the collective velocities and shock stability conditions. The findings are 

confirmed by comparison to Monte Carlo data of a multi-parameter class of simple two-lane driven 
diffusion models, which have the stationary state of a product form on a ring. Going to the hydro- 
dynamic limit, we point out the analogy of our results to the ones known in the theory of differential 
equations of conservation laws. We discuss the singularity problem and find a dissipative term that 
qq . selects the physical solution. 
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! I. INTRODUCTION 

a 

q ■ Driven many-particle systems have been the topic of numerous studies in recent years [1]. Despite 
l— ~~ 1 relatively simple formulation, they have rich dynamic features and phase behaviour and proved to be 
useful testing ground in nonequilibrium physics. Among the new phenomena they highlight are, for 
instance, nonequilibrium boundary- driven phase transitions [2], spontaneous symmetry breaking [3], 
and others. From the mathematical viewpoint, dynamics of many nonequilibrium particle models are 
. Markov processes. The hydrodynamic limit of the latter contributes to the theory of the differential 
^~~] ! equations of conservation laws [4,5]. 

(N| ; If a driven system consists of particles of only one type (one species case), its dynamics can be 
well understood in terms of elementary excitations [6]. Pursuing further the approach of [6], one can 
explain and subsequently predict the stationary phase diagram for systems with arbitrary current- 
density relation [7]. However, multi-species models (i.e. those having two or more different particle 
types, each conserved separately) have so far eluded careful examination. For very recent work, see 
[8]. Earlier studies indicate a phase diagram much richer compared to the one species case, including 
O the existence of new phases and phase transitions [3,9]. The dynamic properties that lead to these 
phase transitions were not studied. 
The present paper presents a step towards understanding of how elementary excitations behave in 
^ ■ a driven system with two species. For that purpose, we propose a new class of models for which the 
stationary state has a simple form on a ring. We study them by analytical means and supplement 
our findings with numerical Monte Carlo simulations and mean field calculations. 

The paper is organized as follows. In section II we define the model and describe the stationary 
state on a ring. In section III we describe how a local perturbation spreads and obtain the eigenvalue 
equation for the collective velocities. In section IV the dynamical evolution of a system with step-like 
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initial condition ( Riemann problem) is considered. We discuss shocks and rarefaction waves and 
their stability. A comparison with partial differential equations obtained by taking the continuum 
limit of the microscopic model is given in section V. Some technical details concerning the derivation 
of the stationary flux and the spreading of a local perturbation can be found in the Appendices. 

II. THE MODEL 

Our model consists of two parallel chains, chain A and chain B. Each chain contains hard-core 
particles which hop randomly to the nearest right or left site if it is empty. Hopping between the 
chains is not allowed as well as occupation of a site by more than one particle (exclusion principle). 
The hopping rate in one chain depends on the configuration of the neighbouring sites in the other 
chain so that for each chain one has the eight possibilities shown in Fig. 1. For the sake of simplicity, 
we consider either the symmetric system (when the hopping rates on the both chains are the same), 
or the antisymmetric system (when the rates on the different chains are left-right reflected). 

The stationary state for a periodic system has a simple form (3) if and only if the hopping rates 
(Fig. 1) satisfy the condition 

a - a L + (3 - (3 L = 2 (7 - 7 l e") = 2(e e" - e L ) (1) 

for the symmetric model and 

a - a L = (3 - (3 L = 7 e~ u - ^ L = e - e L e~ u (2) 

for the antisymmetric one. Here v is an arbitrary real number. For a ring of L sites, the stationary 
probability of the configuration ^1',™^..'.?™^ ^ nen nas ^ ne P rocm ct form 

PnXXn? L = Z ~' II ex P n k m k ). (3) 
fc=i 

Here n^, are occupation numbers for the A- and the 5-chain respectively, i.e., = (n^ = 1), if 
the fc-th site on the A-chain is empty (occupied by a particle). Z is a normalization factor, analogous 
to the partition function in statistical mechanics. One can check that Eqs.(l)-(3) satisfy the station- 
arity requirements by considering gain and loss processes from and to an arbitrary configuration like 
it is done in [10]. One sees from Eq. (3) that different sites k ^ j are uncorrelated. If in addition 
v — 0, then the adjacent pairs of sites are also uncorrelated. 

For demonstration purposes, we shall take general symmetric model (1) and simplify it further: 
firstly, we forbid all backward hoppings and secondly, we set v — 0. With the above restrictions, the 
hopping rates satisfy a + (3 = 2 7 = 2e and can be parametrized by only one parameter (3: 

a = l; 7 = e=(l+/3)/2; a L = (5 L = lL = e L = 0. (4) 

Since the rates must not be negative, (3 is in the range < (3 < 00. (3=1 corresponds to the model 
with no interaction between the chains, called totally asymmetric exclusion process (TASEP) [11,12]. 
The other choice of the rates a = = 7 ^ e, ai = Pl = 1l = £l = was considered in [9]. With 
the choice (4), stationary fluxes j A ,j B of particles on the chains A, B are readily computed to be 

J A (p A ,P B )=P A (1-P A )(1 + (/3-1)p B ) 

j B (p A ,p B )=p B (l-p B )(l + (/5-l)p A ) (5) 

For completeness, we list here exact analytic expressions for the flux in general case (see A for 
details). For symmetric hopping rates (1), 
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J s A ym (p A ,P B ) = (c - lD{e v - 1)/ab/ba + (a - a L - 7 + - P A )/ab (6) 
+ (/?-&- 7 + e L )p A /ba + (7 - e L )p A (l - p A ) 

JsymlP >P ) = JsymlP >P ) 

where /ab (/ba) is the stationary probability to find a particle on chain A (B) and a hole on the 
adjacent site on the other chain. This probability can be obtained from the stationary distribution 

(3), 

Jab 1 (7) 

1 + ^ab — ^ba + y 1 + (^ab — ^ba) — 2JF AB — 2JF BA 

T ST = (l-e-»)p s {l-p T ). 

/ba is obtained by exchange pa Pb in (7). For antisymmetric rates (2) 

Ja A ym (p A ,P B ) = (7 - € L )(e- v - 1)/ab/ba + (a-a L - 1 + e L )(l - p A )f AB (8) 
+ W ~ Pl ~ 7 + ^l)p A /ba + (7 - 6l)p A (1 - P A ) 

Ja B s ym (P A 'P B ) = -^a A s ym (p B >P A ) ) 

/ab being again given by (7). 

III. SPREADING OF A LOCAL PERTURBATION 

We shall study how an initial point-like excitation on the otherwise homogeneous stationary back- 
ground propagates in the system. For our purposes it is convenient to adopt the quantum Hamiltonian 
formulation of the stochastic process. The method is reviewed in detail in [1]. In this formulation, 
the state of our classic stochastic system of two chains of L sites is represented by a vector in a vector 
space, 

i*> c (c*y L ® (c*y L . (9) 

The stochastic dynamics is governed by a quantum Hamiltonian H acting in that vector space, 
d\$)/dt = -H\V), with formal solution \V(t)) = e" m |^(0)). Stationary (i.e., time independent) 
states satisfy evidently H\fy sta , t ) = 0. A particle on site k is represented by the vector (^j and 

a vacancy by the vector at the relevant place in the tensor product (9). E.g., a system with 
only one particle on site k = 2 on the A-chain and no particles on 5-chain is represented as \^/) = 

(0) ® (1) ® (0) 2 ® (0) L ■> e ^ c " Given the state of the system |\&), the expectation value of particle 
density at site k on A-chain is computed by the analogue of the quantum-mechanical formula 

(n k ) = (s\n k \V) (10) 

where (s| = (111 ... 1) is the row vector with all components 1, and h k is a local occupation number 

operator n — ^ ^ ^ ^ acting nontrivially at the k-ih subspace of the tensor product: 

n k = ®n®/ 0i - fc ) J=(jj)- (11) 

Analogously, the expectation value of particle density at site k on 5-chain is computed by averaging 
operator rh k = I® L <g> (J®*-i ® ft ® I® L ~ k ). 
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For our choice of the rates (4) there are no bulk correlations, all configurations with a fixed number 
of particles occur with the same probability (see (3), v — 0). Given the average particle density p z 
on chain Z, the corresponding stationary state within the quantum Hamiltonian formalism is written 
as a product measure 



meaning that there is a probability to find a particle on the A-chain (hk) = p A and a hole (/ — hk) = 
1 — p A at any site k, and (rhk) = p B , (I — rhk) = 1 — p B for the 5-chain. 

We shall study the time evolution of the above homogeneous state in an infinite chain — oo < k < oo 
perturbed at a single site k = 0: 

where \p A p B ) is a stationary state (12). $ A ,$ B are constants, determining the strength and sign of 
perturbation at site 0. We shall see below that only the ratio $ A /$ B is important, hence we consider 
$ A , $ B to be sufficiently small for the averages (14) to be in a physical domain < (no), (rh ) < 1. 
The density profile corresponding to this initial state is given by the average occupation numbers 
(fik) = P A , (rhk) — P B f° r an sites k ^ 0. At the site k — the densities correspondingly are 

(no) = P A + $ A ; (m ) = p B + $ B (14) 

The time evolution of the initial state (13) is given by a Hamiltonian H [1] of the stochastic process 

i*w> = e-«>mo)} = i/p b > + e-«> (* A ^$y + * B ^i^y) i" A " B >- < 15 > 

where we used the stationarity of \p A p B ): H\p A p B ) = 0. 
Consider the quantity 

o/.x (EkK^k- P A )) 1 a u / 1r n 

S " >= (E t (>u--P*)) =iA<E » (16) 

which tracks the position of the center of mass of the excitation on A-chain. Denoting by S the 
change of the center of mass position during the infinitesimal time interval r, we have 

$ A 5 = $ A (S(t + t) - S(t)) « r| £ *<**) = ' E *<5ti " i A ) (17) 

k k 

where we used the lattice continuity equation j^n k = j A ^ 1 — j A . Using the fact that far from the 
excitation there is an unperturbed state with a stationary flux (j A ) = j A , one can shift the summation 
variable in (17) as 

E HjIi - j£) = E(* + i)& A > - E Hjt) = E6' A - A (is) 

A; k k k 

Finally, it can be shown (see an B) that 

so, A (i£» 
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From (17) and (19) we have 

where va = f is the collective velocity of the excitation on the A-chain. Repeating the calculations 
for the £>-chain, we obtain: 

t£* A + w** = " B * B (21) 

Now, if there is an interaction between the chains, the collective velocities must coincide, since the 
perturbation in one chain causes the response in the other and vice versa. Thus vb — va — v, and 
one recognizes in (20,21) the eigenvalue equation X>|$) = v\&), where |$) = (* B ), and V is the 
Jacobian T)^ = dji/dpk- 

The solutions of the eigenvalue problem v™ 11 , fjf 11 and the corresponding eigenfunctions ($2) 
have a transparent physical meaning. Namely, the center of mass of the initial perturbation $ A , $ B 
in the adjacent pair of sites will move with the velocity i£ o11 . An arbitrary initial perturbation |<3>) 
will propagate along the two characteristics vl° n t, v 2° u t. The conserved masses M. z = ^ k (n% — p z ) 
of the splitted components will relate like 0:1/0:2 where 0:1,0:2 are expansion coefficients given by 
|$) = ai|$i) + a 2 \$> 2 )- 

Let us demonstrate the theory in the case of our stochastic model. The Jacobian V ik = dji/dpk is 
readily obtained from (5) 

/(l-2p A )(l + (/?-l)p B ) 03-l) P A (l-p A ) \ 

V (/3-l)P B (l-P B ) (l-2p B )(l + (^-l)p A )J- ^ 

The collective velocities vl° n > v™ 11 are the eigenvalues of V. If f3 — 0, the eigenvalues and corre- 
sponding eigenvectors are given by: 

„,coll 



.coll 



= (l_ p A )(1 _ p B );$i= / 1 V p = Q . (23) 

2 pA_ 2p B + 3p A p B. $ 2 =(J B ) ; (3 = 0. (24) 



Fig. 2 shows the time evolution of the initially perturbed state (13) from Monte Carlo calculations, 
using random sequential update. The background densities are p A = p B = 0.5, which corresponds 
to vl° n = — t>2° u = 0.25, |$i) = (j^), I $2) — (i)- Hence the initial asymmetric excitation (14) 
with <3> A = — $ B must spread to the right, and the symmetric excitation <£> A = <£> B to the left with 
the collective velocities v J oU and -uSf 11 respectively. This is precisely what is seen on the Fig. 2. An 
arbitrary excitation 

' (25) 





$ B / 2 VI/ 2 



will split in two with the masses ratio ($ A — $ B )/($ A + $ B ), spreading apart to the right and to the 
left from the origin. 



IV. SHOCK WAVES, RAREFACTION WAVES AND THEIR COMBINATIONS 

Now we ask the question: if we have prepared the system in a step function (shock) state with 
constant stationary backgrounds (p A , p B ) and (p^, p B ) at the left half space k < and the right half 
space k > respectively, what will happen with their interface? 
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Suppose interface will start moving. Due to mass conservation the Z-component of the interface 
should move with the velocity 

V z (L,R) = 4 Zj % (26) 
Pr ~ Pl 

where we used shortened notation jf^ = j Z (Pl(r)^ Pl(r))i an d the V Z (L,R) marks the fact that 
the velocity is computed between the backgrounds "L" and "R". If V A (L, R) = V B (L,R), the two 
interfaces evolve coherently, similar to the case discussed below (Fig. 4). If however V A (L, R) ^ 
V B (L,R), then the incoherent motion in the A-chain will influence the 5-chain and vice versa, 
destroying the interface. The possible way out for the system is to develop a plateau "0" in the 
middle, interpolating between the plateaus "L" and "R", as shown on Fig. 4. Consequently, instead 
of one there will be two interfaces in each chain: the interface L\0 between "L" and "0" and the 
interface 0|R. We must require the velocities in the A- and I?-chains to be the same. The interface 
L\0 has the velocity 

y(L ' 0) - p a l -pz ~ ft-fi (27) 

and analogously for the interface 0|R 

vm m - tW^tM - j B ( fi )-J B (°) , 2R , 
m R) ~ ~ p?-p„ b (28) 

The solutions of Eqs. (27,28) define the location of possible middle plateau densities p^,p B . Since 
(27,28) are nonlinear, they can have several solutions or no solutions at all. If a solution exists (see 
Fig. 3), one must require additionally V(L, 0) < V(0, R) because the plateau "0" has to expand, and 
check the shock stability as discussed below. 

In order to study the shock stability let us consider a shock of the form u L\0\R" consisting of three 
consecutive plateaus at densities pf , p% and A small deviation at the plateau U K" (K = 0, L, R) 
will split into two local excitations with the velocities v f n (K) > V2° n (K) as discussed in section III. 
In order for the shock to be stable, all local excitations have to be absorbed by the interfaces, that 
is, 



■ coll (L), vf\L) > V(L, 0) > vf\0) (29) 



for the interface L|0 and 



< u (0) > V(0, R) > v c 2 ol \R),vf\R) (30) 



for the interface 0|R. An example of such a double shock is shown on Fig. 4. The densities of particles 
PoSPo on t ne middle plateau satisfy Eqs. (27,28), which are graphically solved on Fig. 3. There are 
two solutions, one of which is realized (Fig. 4), while the other one violates V(L, 0) < V(0,R). One 
can check that (29,30) are satisfied. 
If the second shock condition (30) is not satisfied, but instead one has (30) 

< u (0) < V(0,R) < vf\R),vf\R); vf\0) > V(0,R), (31) 

then the local perturbations will destroy the sharp interface, leading to a rarefaction wave connecting 
the plateaus 0|R, similar as discussed for one-component systems [1,7]. One may ask what happens 
if neither shock-wave nor rarefaction-wave condition are satisfied. In this case a combination of both 
shock and rarefaction wave may be formed as illustrated on the example Fig. 5. There for simplicity 
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we have taken symmetric initial conditions so that Eqs (27,28) are satisfied for arbitrary p A = p B - The 
velocity of the mass transfer (26) between the left (L) and right (R) plateaux respectively is the same 
for both chains V A (L,R) = V B (L,R), but the collective velocities v™ l \R), v^ ol \R) < V(L,R) and 
vl° u (L), vf ll (L) < V(L,R) satisfy neither of shock-type (30) nor of rarefaction-wave (31) criterium. 
As a result, the compromise is made: part of the profile with the high densities p C (p L , p*) develops 
a shock while for the small densities p C (p*, Pr) the rarefaction wave is formed, see Fig. 5. Indeed for 
the interface {pl,P* + e ), e ^ 1, the shock condition analogous to (30) is satisfied while the interface 
(p* — e, pr) satisfies respective rarefaction wave condition (31). The level p* = 0.525 is defined by 
the crossing point V(pl,P*) = vf ll (p*) and can be predicted by hypothetical consideration of the 
initial condition as a sequence of small plateaus (shocks) at each level of density. All small shocks 
above p = p* condense in a single shock while those below p* form the rarefaction wave. Similar 
analysis can be performed for other initial conditions. Note that for certain class of initial conditions 
we observe shocks of even more complicated structure. Their analysis will be presented elsewhere. 



V. HYDRODYNAMIC LIMIT: COMPARISON WITH THE THEORY OF PARTIAL DIFFERENTIAL 

EQUATIONS 



For most notions we have discussed in the framework of the stochastic particle system, one can find 
the respective analogies in the theory of partial differential equations (PDE). The naive continuum 
(Eulerian) limit of our stochastic dynamics on the lattice hk(t) — > p A (x,t), rh k {t) — > p B (x,t) is a 
system of conservation laws 

where j z is given by Eqs(5-8). Here and below in this section we shall use p z (x, t) for a continuously 
changing variable, not to be confused with constant p A ,p B from section III. Such systems of 
conservation laws are studied e.g. in [13,14]. 

The eigenvalues of the Jacobian |^ (playing the role of collective velocities) are the characteristic 
velocities. For scalar conservation law dp/dt + dj(p)/dx = with dj(p)/dp = v(p) it follows that 
the line x = v(p)t, called characteristic, defines the space-time trajectory on which the local density 
p(x,t) stays constant. For systems (32) of conservation laws, the situation is more complicated. 
However also there one can find two functions u>i(p A , p B ); i = 1, 2, called Riemann invariants, which 
are constant along the respective characteristics ^ = Vi [14]. 

Consider now a shock of the type drawn in Fig. 4. In the hydrodynamic limit, the interface 
region between the plateaus of constant densities will squeeze to a single point, giving rise to dis- 
continuous change. Discontinuities in PDE theory are known to satisfy the so-called jump condition 
v s (p+ — p z ) = (j+ — j z ) where F + and F_ are the values of function F in the right and left edges 
of the discontinuity, and v s is the speed of the propagation of the discontinuity. Comparing with 
(26) we recognize in v s the shock velocity. 

It is well known that an arbitrarily chosen smooth initial profile p z (x, 0) will develop a singularity 
after finite time t [13,14]. To cure the singularity problem for the PDE, the simplest possible approach 
suggests adding a vanishing viscosity term (k^^-; k — > o) to the right-hand side of (32). This is 
enough to avoid singularities and by numerical integration we find that this regularization term leads 
to the correct answer for the initial Riemann problem, as compared to the stochastic model. 

Another possibility to obtain a viscosity term is to average exact lattice continuity equations of 
the stochastic process 

^k = jLx-k (33) 
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^rn k =jti-jf (34) 

for occupation number operators (ri k ) — > p A (x,t), (m k ) — > p B (x,t), allowing for continuous change of 
space k, k + 1 — » x, x + dx. For the case (4), the flux operator j£ can be obtained from the general 
expression (Al) and it reads 

3k = M 1 ~ n k +i) (l + (™>k + m k+1 )j . (35) 

j B is obtained by an exchange n <-> rh in the above. We substitute (35) into (33), (34), average, 
factorize and Taylor expand the latter with respect to site spacing dx as e.g. (rh k+ i) = p B (x,t) + 
dx 9p g ^'^ + (dx 2 /2) d gjf'^ + Keeping the terms up to dx 2 in the expansion, we obtain 

dx d d . . 

K = T 0; « = 2K *' (38) 

where j z (p A ,p B ) are given by the (5). 

We found by a numerical integration of (36), (37) that also here the correct result are obtained for 
the step-function initial conditions. It seems therefore that the choice of viscosity matrix is rather 
arbitrary, if initial step-function conditions are chosen. However, the PDE becomes more sensitive 
if solved on a finite interval with fixed boundary values. In this setting, the choice of viscosity is 
important, since different choices give different answers. The details will be published elsewhere. 



VI. CONCLUSION 



To conclude, we have studied a two lane particle exclusion process on the microscopic level, focusing 
on the temporal behaviour of the elementary local excitations, domain walls (shocks) and rarefaction 
waves. By analyzing the flow of localized excitations and calculating their collective velocities we 
derived a criterion for the stability of shocks, somewhat analogous to our recent analysis of systems 
with a single conservation law [1]. However, unlike in systems with a single conservation law, shocks 
generically come in pairs, since the two conserved densities give rise to two distinct collective veloc- 
ities. Because of the interaction between the chains, or more generally, between the two conserved 
densities, these velocities are the eigenvalues of the Jacobian of the current- density relation. The 
eigenvectors of the Jacobian parametrize the expansion of the strength of a generic excitation into 
the two eigenmodes of the systems. These eigenmodes (corresponding to the special case of a single 
excitation) correspond to a special tuning of the strength of the excitation in each conserved density: 
The relative strength for an eigenmode with fixed collective velocity is the ratio of the components 
of the corresponding eigenvector. Initial profile not satisfying the stability criterion for shocks evolve 
into rarefaction waves or more complicated structures. Since nowhere in our analysis we make use 
of the specific properties of our model we argue that as in systems with a single conservation law, 
all the properties discussed above can be derived from the macroscopic current, irrespective of the 
microscopic details of the model. 

Going one step further we take a naive continuum limit (Euler scale) to obtain a system of coupled 
nonlinear PDE's. Thus we obtain a microscopic interpretation for the characteristics (as flow of 
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localized perturbations) and for the jump condition for shock solutions. Monte-Carlo simulation of 
the model as well as numerical integration of the PDE's suggest that the uniqueness problem for 
the Riemann problem can be resolved by using a quite arbitrary viscosity matrix with vanishing 
viscosity. However, for the stationary solution with fixed boundary values the problem appears 
to be more intricate. A detailed analysis is necessary and will be presented in future work. The 
hydrodynamic limit of another family of lattice gas models with two conservation laws, differing 
from ours by internal symmetries, has been studied recently [8]. They give rise to a different set of 
PDE's, but we believe that our analysis can be applied to this family as well. On the other hand we 
also expect that the mathematically rigorous work of [8] can be generalized to models of the type 
considered here. 
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APPENDIX A: CURRENT-DENSITY RELATION IN THE GENERAL CASE 

By definition, the stationary flux is written as sum of the hopping rates times the stationary 
probabilities n of the corresponding local configurations, see Fig. 1: 

j A ( P A , P B ) = an Q + pn CD + 7 n a + en a 

- a L n a - (3 L n (::) - lL n a - e L n a (Ai) 

Here we introduced the short notation n (°°) for the probability to find the configuration with one 
particle (filled circle •) and 3 holes (empty circle o), arranged like in Fig. 1, first configuration on the 
upper row), in a steady state with average densities p A and p B . Analogously, n (") is the probability 
to find 3 particles, 1 hole as in Fig. 1, second configuration on upper row. In terms of occupation 
number operators ri k ,m k for chains A (bottom circles) and B (upper circles), 

^ (") = ("fc(l - n k+ i)m k m k+ i) , (A2) 

and so on. If the rates satisfy (1) and (2) in the symmetric and antisymmetric case respectively, then 
the correlation function can be factorized due to (3) e.g. 

(n k (l - n k+1 )m k m k+1 ) = (n k m k )(((l - h k+1 )m k+1 ) = p B (h k m k ) - (h k m k ) 2 . (A3) 

Above we used the translational invariance and the fact that (rh k ) = p B , (n k ) = p A . Factorizing 
Eq.(Al), and using (1), (2), one obtains 

j a (p a ,p b ) = k none) + (i-p A )(a-a L - 1 + e L )no 

+ p A ((3 - f3 L - 7 + cl) n C) + P A (1 - P A )(7 - cl) 

where 

r, , r> f (e — 7L)(e 1/ — 1), symmetric case , . 

K = -- + a L -(3 + (3 L + e + 1 -e L - lL = ^ _]^ u _ ^ case (A4) 

Finally, n (°) = (n k {l — m k )) can be calculated directly from the stationary distribution (3) which 
after some algebra gives the expression (7). n (•) is obtained from (7) by exchanging p A <-> p B . 
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APPENDIX B: PROOF OF THE EQ.(ll) 

Consider C = J2k(j A ~ i A )> which is by definition 



£ = E<*l(i A -J A K 



Ht 



A n ~ P 



B m 0- P 



pA(l_ p A ) pB (1 _ pB) ^ 



(Bl) 



(s\ is the constant row vector (1 1 1 ... 1) (see [1] for details). First, (s\j A \p A p B ) = j A by the definition 
of the stationary flux. Thus (Bl) is simplified as 



:-a -tft^A n o - P 



pA(l-pA) 



^A-Ht^B m P I A B 



pB(l-pB) 



Due to translational invariance, the above expression can be rewritten as 



V ° p a (i-p a ) ! 



P A P B ) + (slj^e 



B*,B„-Ht m k P \p A p B } 



P B (1-P B ) 



(B2) 



(B3) 



Because the total number of particles in each chain is conserved, the Hamiltonian H commutes with 
J^k^k, J2k™k- Using this, and the fact that the |p A p B ) is stationary, the term e~ Ht can be deleted 
from (B3). Substituting the definition of |p A p B ) from (12) into the expression below, we have 

E(^-p a )IpV) 



E 

k 



v P A 



l-p* 



i-P 

, p a 



i-P 

, P B 

d 



B x 



A „B\ 



„A/-i „A\ " I A „ 

p 0--P )^a> p 



Here we have used the explicit representation of the particle number operator h 




1 



(B4) 



Anal- 



ogously, Efc (/ho - P B ) \P A P B ) = P B (1 - P B )^b \P A P B }- Substituting this together with (B4 )in (B3) 
we obtain Eq.(19). 

Note the specific fact that the stationary state in our system is a product measure (12). However 
the validity of the result (19) extends to a much wider class of driven systems with short-ranged 
correlations. 
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FIG. 1. The eight allowed elementary hopping processes for the first chain, and their rates. In the study we choose the rates 
for the second chain either to be the same, Eq.(l) or to be antisymmetric (reflected with respect to the first chain), see Eq.(2). 
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FIG. 2. Time evolution of a point-like initial perturbation (+5, ±6) on a 2-chain driven system. The parameters are: the 
background densities p A = p B = 0.5, /? = 0. LEFT graph: at t = 0, the asymmetric perturbation Sp A = —Sp B = 0.5 is put at 
the middle site 150. The average density profiles after t = 100 and t = 200 Monte Carlo evolution steps are shown, averaged 
over 6 * 10 different histories. The component A is depicted with points and the component B with lines. RIGHT graph: initial 
perturbation is symmetric 8p A = 8p B = 0.5. B- component evolution (not shown) is identical to the one of the A-component. 
The asymmetric perturbation moves with collective velocity 0.25 to the right, and the symmetric one to the left, in accordance 
with the theory (see section III) . 
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3. Curves, showing locus of the points p A ,p B , solving Eq.(27) with p A 



0.6, pf = 0.75 (bold curves). (3 = 0.2. 



0.15, pi = 0.3 (thin curves) and Eq.(28) with 
Crossings of the bold curves with the thin curves indicate possible solutions of 



(27,28). 
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t = 200 
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FIG. 4. Formation of a shock wave. Parameters are: j3 = 0.2. LEFT: the initial distribution: ^4-particles are distributed 
randomly with an average densities 0.15(0.6) on the left (on the right). Corresponding B-particles initial densities are 0.3(0.75). 
RIGHT: Result of Monte Carlo evolution after t = 200 Monte Carlo steps, averaged over 2 * 10 5 different histories. The density 
profiles of A and B particles are depicted with points of different sizes, and the lines with theoretically expected middle shock 
values p A « 0.3146, p B » 0.8323 are drawn for comparison. 
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FIG. 5. Time evolution of a step-like initial profile, leading to shock wave and rarefaction wave coexistence. The parameters 
are: the left densities p A = p B = 0.95, the densities on the right p A = p B = 0.3, (3 — 0. The graph shows the average 
density profiles at t = 0, 100, 150, 200, symmetric in both components. Average over 10° histories is made. For the densities 
p = p* = 0.525 and higher, the shock wave is formed, while for the lower densities p < p* one observes the rarefaction wave 
(see the end of section IV for details) . The level p = p* is marked by a thin line. 
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